home *** CD-ROM | disk | FTP | other *** search
/ Games of Daze / Infomagic - Games of Daze (Summer 1995) (Disc 1 of 2).iso / x2ftp / msdos / math / random-c / random.c < prev   
C/C++ Source or Header  |  1990-05-28  |  5KB  |  163 lines
  1. /*
  2. C This random number generator originally appeared in "Toward a Universal 
  3. C Random Number Generator" by George Marsaglia and Arif Zaman. 
  4. C Florida State University Report: FSU-SCRI-87-50 (1987)
  5. C 
  6. C It was later modified by F. James and published in "A Review of Pseudo-
  7. C random Number Generators" 
  8. C 
  9. C THIS IS THE BEST KNOWN RANDOM NUMBER GENERATOR AVAILABLE.
  10. C       (However, a newly discovered technique can yield 
  11. C         a period of 10^600. But that is still in the development stage.)
  12. C
  13. C It passes ALL of the tests for random number generators and has a period 
  14. C   of 2^144, is completely portable (gives bit identical results on all 
  15. C   machines with at least 24-bit mantissas in the floating point 
  16. C   representation). 
  17. C 
  18. C The algorithm is a combination of a Fibonacci sequence (with lags of 97
  19. C   and 33, and operation "subtraction plus one, modulo one") and an 
  20. C   "arithmetic sequence" (using subtraction).
  21. C======================================================================== 
  22. This C language version was written by Jim Butler, and was based on a
  23. FORTRAN program posted by David LaSalle of Florida State University.
  24. */
  25.  
  26. #include <stdio.h>
  27. #include <stdlib.h>
  28.  
  29. #define TRUE -1
  30. #define FALSE 0
  31. #define boolean int
  32.  
  33. static void rmarin(int ij, int kl);
  34. void ranmar(float rvec[], int len);
  35.  
  36. void main()
  37. {
  38.     float temp[101];
  39.     int ij, kl, len, i;
  40.  
  41.     /* random seeds for the test case: */
  42.     ij = 1802;
  43.     kl = 9373;
  44.  
  45.     /* do the initialization */
  46.     rmarin(ij,kl);
  47.     
  48.     /* generate 20,000 random numbers */
  49.     len = 100;
  50.     for (i=1; i<=200; i++)
  51.         ranmar(temp, len);
  52.         
  53. /*
  54. C If the random number generator is working properly, the next six random
  55. C numbers should be:
  56. C           6533892.0  14220222.0  7275067.0
  57. C           6172232.0  8354498.0   10633180.0
  58. */
  59.      len = 6;
  60.      ranmar(temp,len);
  61.      
  62.      for (i=1; i<=6; i++)
  63.          printf("%12.1f ", 4096.0*4096.0*temp[i]);
  64. }
  65.  
  66. static float u[98], c, cd, cm;
  67. static int i97, j97;
  68. static boolean test = FALSE;
  69.  
  70. static void rmarin(ij,kl)
  71. int ij, kl;
  72. {
  73. /*
  74. C This is the initialization routine for the random number generator RANMAR()
  75. C NOTE: The seed variables can have values between:    0 <= IJ <= 31328
  76. C                                                      0 <= KL <= 30081
  77. C The random number sequences created by these two seeds are of sufficient 
  78. C length to complete an entire calculation with. For example, if sveral 
  79. C different groups are working on different parts of the same calculation,
  80. C each group could be assigned its own IJ seed. This would leave each group
  81. C with 30000 choices for the second seed. That is to say, this random 
  82. C number generator can create 900 million different subsequences -- with 
  83. C each subsequence having a length of approximately 10^30.
  84. C 
  85. C Use IJ = 1802 & KL = 9373 to test the random number generator. The
  86. C subroutine RANMAR should be used to generate 20000 random numbers.
  87. C Then display the next six random numbers generated multiplied by 4096*4096
  88. C If the random number generator is working properly, the random numbers
  89. C should be:
  90. C           6533892.0  14220222.0  7275067.0
  91. C           6172232.0  8354498.0   10633180.0
  92. */
  93.     int i, j, k, l, ii, jj, m;
  94.     float s, t;
  95.     
  96.     if (ij<0 || ij>31328 || kl<0 || kl>30081) {
  97.         puts("The first random number seed must have a value between 0 and 31328.");
  98.         puts("The second seed must have a value between 0 and 30081.");
  99.         exit(1);
  100.     }
  101.     
  102.     i = (ij/177)%177 + 2;
  103.     j = ij%177 + 2;
  104.     k = (kl/169)%178 + 1;
  105.     l = kl%169;
  106.     
  107.     for (ii=1; ii<=97; ii++) {
  108.         s = 0.0;
  109.         t = 0.5;
  110.         for (jj=1; jj<=24; jj++) {
  111.             m = (((i*j)%179)*k) % 179;
  112.             i = j;
  113.             j = k;
  114.             k = m;
  115.             l = (53*l + 1) % 169;
  116.             if ((l*m)%64 >= 32) s += t;
  117.             t *= 0.5;
  118.         }
  119.         u[ii] = s;
  120.     }
  121.     
  122.     c = 362436.0 / 16777216.0;
  123.     cd = 7654321.0 / 16777216.0;
  124.     cm = 16777213.0 / 16777216.0;
  125.     
  126.     i97 = 97;
  127.     j97 = 33;
  128.     
  129.     test = TRUE;
  130. }
  131.  
  132. void ranmar(rvec, len)
  133. float rvec[];    /* len random numbers are placed in rvec[1..len] */
  134. int len;
  135. /*
  136. C This is the random number generator proposed by George Marsaglia in 
  137. C Florida State University Report: FSU-SCRI-87-50
  138. C It was slightly modified by F. James to produce an array of pseudorandom
  139. C numbers.
  140. */
  141. {
  142.     int ivec;
  143.     float uni;
  144.     
  145.     if (test==FALSE) {
  146.         puts("Call the init routine rmarin() before calling ranmar().");
  147.         exit(2);
  148.     }
  149.     for (ivec=1; ivec<=len; ivec++) {
  150.         uni = u[i97] - u[j97];
  151.         if (uni < 0.0) uni += 1.0;
  152.         u[i97] = uni;
  153.         i97--;
  154.         if (i97==0) i97 = 97;
  155.         j97--;
  156.         if (j97==0) j97 = 97;
  157.         c -= cd;
  158.         if (c<0.0) c += cm;
  159.         uni -= c;
  160.         if (uni<0.0) uni += 1.0;
  161.         rvec[ivec] = uni;
  162.     }
  163. }